clear all;
close all;
delete(gcp('nocreate'));
% Current folder for Matlab codes
fp.matlab = pwd;
common_path = extractBefore(fp.matlab,"\matlab_dir");

% Paper folder;
fp.paper = sprintf('%s\\paper', common_path);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We estimate the model at lambda=0.95;
fp.lambda = 0.95;

%model primitives;
fp.ages =9:12;
fp.periods=length(fp.ages);
fp.percentiles_kid = linspace(5,95,10);
fp.hc_points_kid = length(fp.percentiles_kid);

fp.Min_Time = 0.01;
fp.Max_Time= 1;
fp.inv_points=20;
fp.percentiles_inv= linspace(5,95,fp.inv_points);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%indexes for latent data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.ind_data.child_id = 1;
fp.ind_data.child_age = 2;
fp.ind_data.n_sim = 3;
fp.ind_data.child_C = 4;
fp.ind_data.child_C_tp1 = 5;
fp.ind_data.inv = 6;
fp.ind_data.peers = 7;
fp.ind_data.peers_log = 8;
fp.ind_data.ParStyle = 9;
fp.ind_data.peers_tp1 = 10;
fp.ind_data.peers_log_tp1 = 11;
fp.ind_data.neighborhood = 12;
fp.ind_data.moved = 13;

fp.ind_data_network.child_id=1;
fp.ind_data_network.child_age=2;
fp.ind_data_network.n_sim=3;
fp.ind_data_network.child_C=4;
fp.ind_data_network.n_friends=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loading Estimated Moments;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

initial_conditions = importdata('initial_conditions.txt');    
fp.schools_distrib = initial_conditions;

mean_income_by_neighborhoods = importdata('family_income_neighborhoods.txt');
mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
fp.mean_income_by_neighborhoods = mean_income_by_neighborhoods;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%random mumber draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rng(10); %fix seed;
fp.n_rip_sim = 200;                         % # of simulation of the economy;
fp.tot_draws_mcintegration_network = 50 ;   % # of simulation of the network;

fp.init_draws{1} = norminv(rand(fp.schools_distrib(1,3),1,fp.n_rip_sim));
fp.init_draws{2} = norminv(rand(fp.schools_distrib(2,3),1,fp.n_rip_sim));
fp.init_draws{3} = norminv(rand(fp.schools_distrib(3,3),1,fp.n_rip_sim));
fp.init_draws{4} = norminv(rand(fp.schools_distrib(4,3)+50,1,fp.n_rip_sim));

fp.peers_backward{1} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(1,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{2} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(2,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{3} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(3,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{4} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(4,3)+50,fp.tot_draws_mcintegration_network, fp.periods-1);

fp.peers_forward{1} = rand(fp.schools_distrib(1,3),fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{2} = rand(fp.schools_distrib(2,3),fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{3} = rand(fp.schools_distrib(3,3),fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);

fp.PS_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COUNTERFACTUALS with general 2 CES case;

%Figure E-1: Policy and Scaling Effects on Skills (Alternative Calibration)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.origin_neigh         = 1;
fp.receiving_neigh      = 4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate Baseline Economy        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('param0_2CES_calibrated');   
%Calibrated Parameters with a second CES technology (for authoritarian
%parents) with the following parameterization (see equation (12) for the general CES specification): 
% alpha_1,1=0.7221, alpha_2,1=0.6734, alpha_3,1=0.30, alpha_4,1=-0.40, alpha_5,1=0.7160. 
% All the rest of other model`s parameters are kept at the estimated values. 
%This new parametrization generates, compared to the Cobb–Douglas case,
%higher substitutability of time vs peers, and higher complementarity with
%current stock of skills; 

fp.do_estimation =0;
fp.do_counter = 0;            
fp.do_counter_type = 0 ;

obj_eval = obj_function( param0_2CES_calibrated, fp );
data_baseline_origin = obj_eval.simul_data{fp.origin_neigh};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Counterfactuals;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.peers_backward_original = fp.peers_backward;

%Here we double the shocks for peers formation in the backward problem
%because in the counterfactuals we double the state space (moving families
%with different preferences into the new neighborhood);
for j = 1:1:4
   
    fp.peers_backward{j} = [ fp.peers_backward{j} ; fp.peers_backward{j} ] ; 
 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Policy #2: Policy Effects of moving average many children from Sending Neighborhood);

%Simulate Counterfactual Policy:
%                               1. Counterfactual Data
%                               2. Policy Effects (Counterfactual VS. Baseline)  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fp.max_kids_moved  = 90;
fp.min_kids_moved  = 41;
fp.do_counter = 1;            
fp.do_counter_type = 2;

%Load baseline results;
load('pe_skills_receiving')
load('pe_skills_moved')

pe_skills_moved_2CES      = NaN(1,fp.max_kids_moved);
pe_PS_moved_2CES          = NaN(1,fp.max_kids_moved);
pe_inv_moved_2CES         = NaN(1,fp.max_kids_moved);
pe_skills_receiving_2CES  = NaN(1,fp.max_kids_moved);
pe_PS_receiving_2CES      = NaN(1,fp.max_kids_moved);
pe_inv_receiving_2CES     = NaN(1,fp.max_kids_moved);

data_baseline = obj_eval.simul_data;


set_children_to_move = NaN(length(fp.min_kids_moved:1:fp.max_kids_moved),fp.n_rip_sim);

for h = 1 : fp.n_rip_sim

set_children_to_move(:,h) = randsample(fp.min_kids_moved:1:fp.max_kids_moved,length(fp.min_kids_moved:1:fp.max_kids_moved));

end

ind = 1;
for n_kids =  [1, 5:5:50]
    
    n_kids
    
    fp.n_kids_moved    = n_kids;

    fp.shock_child_moved = set_children_to_move(1:n_kids,:);

    cd(fp.matlab)
    obj_counter2 = obj_function( param0_2CES_calibrated, fp );
    
    data_counter2 = obj_counter2.simul_data;
    fp.id_moved_children = obj_counter2.id_moved_children;

    pe = policy_effects_counter2(data_counter2,data_baseline,fp);
    
    pe_skills_moved_2CES(ind)=pe.pe_skills_moved;
	pe_PS_moved_2CES(ind) = pe.pe_PS_moved;
	pe_inv_moved_2CES(ind) = pe.pe_inv_moved;
    pe_skills_receiving_2CES(ind) = pe.pe_skills_receiving;
	pe_PS_receiving_2CES(ind) = pe.pe_PS_receiving;
    pe_inv_receiving_2CES(ind) = pe.pe_inv_receiving;

    ind = ind + 1;
           
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Figure E-1: Policy and Scaling Effects on Skills (Alternative Calibration)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

grid_x= [1, 5:5:50];
grid_y = grid_x;
cd(fp.paper)

h1=figure;
scatter(grid_x ,pe_skills_moved(1:length(grid_y)) ,  'b')
hold on
p = polyfit(grid_x ,pe_skills_moved(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
scatter(grid_x ,pe_skills_moved_2CES(1:length(grid_y)) ,  'r')
p = polyfit(grid_x ,pe_skills_moved_2CES(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)    
    
legend('Policy Effects (Baseline)','Quadratic Approx', 'Policy Effects (Alternative Calibration)' , 'Quadratic Approx' )
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Moved Children)');
set(gca,'XTick',0:5:length(fp.min_kids_moved:1:fp.max_kids_moved))
set(gca,'YTick',0.06:0.02:0.24)
set(gca,'YLim', [0.06 0.24 ])
%hgexport(h1, 'Count1_calibrated_2CES.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h1, 'Count1_calibrated_2CES.eps',hgexport('factorystyle'), 'Format', 'eps');

h4=figure;
scatter(grid_x ,pe_skills_receiving(1:length(grid_y)) ,  'b')
hold on
p = polyfit(grid_x ,pe_skills_receiving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
scatter(grid_x ,pe_skills_receiving_2CES(1:length(grid_y)) ,  'r')
p = polyfit(grid_x ,pe_skills_receiving_2CES(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10) 
legend('Policy Effects (Baseline)','Quadratic Approx', 'Policy Effects (Alternative Calibration)' , 'Quadratic Approx' )
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Receiving Children)');
%hgexport(h4, 'Count4_calibrated_2CES.png',hgexport('factorystyle'), 'Format', 'png');    
%hgexport(h4, 'Count4_calibrated_2CES.eps',hgexport('factorystyle'), 'Format', 'eps');     

  
    
cd(fp.matlab)

